SPH Simulations of Accretion Disks and 

Narrow Rings 



Sarah T. Maddison 1 * James R. Miirray 2 ^and Joe J. Monaghan 

1 Maths Department, Monash University, Clayton, 3168, Australia 
2 CITA, University of Toronto, Ontario, M5S 1A7, Canada 

June 6, 1995 



Abstract 

We model a massless viscous disk using Smoothed Particle Hydro- 
dynamics (SPH) and note that it evolves according to the Lynden-Bell 
& Pringle theory (1974) until a non-axisymmetric instability develops 
at the inner edge of the disk. This instability may have the same ori- 
gin as the instability of initially axisymmetric viscous disks discussed 
by Lyubarskij et al. (1994). To clarify the evolution we evolved single 
and double rings of particles. It is actually inconsistent with the SPH 
scheme to set up a single ring as an initial condition because SPH 
assumes a smoothed initial state. As would be expected from an SPH 
simulation, the ring rapidly breaks up into a band. We analyse the 
stability of the ring and show that the predictions are confirmed by 
the simulation. 



1 Introduction 

As a test case for a two dimensional SPH code (for an overview of SPH, 
see for example Monaghan 1992) developed to study the viscous evolution 
of accretion disks we modelled the spreading of a narrow viscously shearing 
ring of matter. We wished to compare simulations with the results of thin 
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disk theory as developed by Lynden-Bell and Pringle (1974). Lubow (1991) 
used the rate of spread of an axisymmetric ring of viscously interacting SPH 
particles to estimate the shear viscosity due to the artificial viscosity term 
in his code. Flebbe et al. (1994) tested a tensor form of a general Navier 
Stokes viscosity term by comparing an axisymmetric ring simulation with 
Lynden-Bell and Pringle's result for the evolution of a 8 function density 
distribution. Because SPH cannot be used with a 8 function density, Flebbe 
et al. used an approximate initial condition. We prefer to be consistent 
and simulate a configuration with an initial state which can be resolved by 
SPH. In particular we looked at an axially symmetric annulus with an initial 
surface density 

E (r) = ( exp{-(^) 2 } n < r < r 2 (1) 
I elsewhere 

where r and / are the radius of maximum density and width of the Gaussian 
respectively. The annulus is assumed to be Keplerian with the only forces 
considered being the gravitational attraction of the central object and bulk 
and shear viscosity forces within the annulus. Pressure and disk self-gravity 
are not considered. To simplify the analysis a constant kinematic viscosity 
v was assumed. For the above initial condition the general solution for the 
surface density at radius r and time t is given by 



£(M) 




where /i/ 4 (x) is a modified Bessel function of the first kind. For large argu- 
ments, Ii/ 4 (x) exp (x) /\/2irx. Thus for small t we have the approximate 
solution: 

{ (r'-r) 2 ] 

Following the equations of Pringle (1981), the simulation is given an initial 
radial velocity 

where E is our initial Gaussian. 
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Figure 1: The time evolution of the surface density of a Gaussian 
disk. The solid lines denote the theoretical solution at the times 
shown. The heavy points show the corresponding SPH results. For 
this simulation we used r\ = 0.80, r$ = 0.85, = 0.90, / = 0.025 
and v = 2.5 X 10" 4 . 



The theoretical and simulated spread of an annulus with a Gaussian density 
profile is shown in figure 1. The best results were obtained when the parti- 
cles were set up in concentric rings so as to give a uniform number density 
throughout the annulus. The Gaussian profile was set up by giving each par- 
ticle a mass proportional to S (r). In this particular case 22,420 particles in 
21 rings were used (see figure caption for further parameters of the model). 

After ten revolutions of the particles at r = the inner edge of the disk 
started to break up. Other simulations of disks made up of more rings gave 
similar results. After several rotation periods, the viscosity spread the rings 
to the point where the innermost ring became somewhat separated from the 
remaining rings. It then started to lose shape and break up. This effect 
marched its way through the remaining disk. One would like to know if 
this is a genuine instability or an artifact of the numerical method. In the 
following we give a partial answer by analysing the stability of a single ring. 
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2 Stability Analysis 



We consider the stability of a single viscous ring of equi-spaced SPH par- 
ticles in orbit with angular velocity fi around a central mass M. Again, 
gravitational and pressure forces of the ring particles have been neglected to 
investigate the viscous forces. 

Each particle i has unperturbed and perturbed positions in radial and az- 
imuthal coordinates: 



where a is the ring radius, fit is a rotating reference line, © 8 is the un- 
perturbed azimuthal position, and Oi and aqi are the azimuthal and radial 
perturbations respectively. 

The radial and azimuthal equations of motion governing a particle in the ring 
are: 



(r, <f>) 



(a,Slt + Qi) 

(a(l + qi),nt + Qi + <Ti) 





d 2 cf> dr dcf> 



— - + 2 = 

dt 2 dt dt 



(4) 



where F = (F r} F^) is the viscous force per unit mass. 




dr 2 H dr QM 



(5) 



d 2 a dq a 2 



dT 2+2 dr QM 



(6) 



For the non- viscous case, where F r = = 0, 



q — 3q — 2(7 
v + 2q 







(7) 



. 




We consider perturbations of the form: 

q = A exp {i(k® + lot)} 
a = B exp {i(k® + lot)} , 



(9) 
(10) 
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where u is the frequency of the disturbance and k is the wave number. Sub- 
stituting (9) and (10) into the linearised equations (7) and (8), we find that 
to = 0, ±1, with a double root at u = 0. The general solution for this non- 
viscous case is the superposition of (a) an oscillation, (b) a shift in © 8 without 
changing a, and (c) the motion associated with all particles being shifted to 
a different circular orbit (see Murray 1994). 

Returning now to the viscous case, we need the SPH viscous forces. Assuming 
that the SPH kernel, Wij, is Gaussian (normalised for two dimensions) we 
can use: 

where h is the SPH smoothing length and r 8J = r 8 - — r r Substituting this into 
the formula for the viscous force per unit mass (Monaghan 1992, eqn 4.1) and 
using pi = errii/h 2 , where pi and ra 4 - are the density and mass respectively 
of particle z, and e is a positive dimensionless constant of order unity that 
depends on the spacing of the particles, 



2^ — acijh 
j Pi 
2ach » , 
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rf; + f] 2 



rijWij, 



where a is the viscosity coefficient, c is the speed of sound (replacing c 8J 
since the sound speed is constant in the isothermal case) and rj prevents 
singularities in the viscous force for very small particle separation. We can 
ignore rj for the purpose of this stability analysis. The kinematic viscosity is 
given by v = ahc/8. 

Since we are considering the linearised equations, we need only consider per- 
turbations to Vij.Tij since it vanishes in the unperturbed case. After lineari- 
sation (see appendix) we find: 

Vij.rij ~ (<Tj - 6-i) sin( 
r% ~ 2(1 - cos C), 

where ( = Qj — The scaled radial component of the force per unit mass 
on particle i is then given by: 

,2 



F r = -^(& j -& i )smCW ij , (11) 



QM 2 
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and the tangential force per unit mass component on particle i is: 



2 

^ = !E(^-^(l + cosC)H4-, (12) 



QM 2 
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where g = 2ach/e. 



We now have the full SPH formulation of the equation of motion of a viscous 
ring, by combining (5) & (6) and (11) & (12), and substituting (9) and (10). 
The summation is over all N particles in the ring. However, only nearby 
particles contribute to the sum, which can thus be replaced by an infinite 
sum. It is then clear that the odd components of these sums must vanish, 
leaving the equation of motion as: 

(-u 2 -3)A-2iuB = gBu^ (13) 

-u 2 B + 2iu)A = -gBiu$ , (14) 

where 

* = ^£ s m(K)shiCH4- 

3 

$ = i]T(l -cos(K))(l + COS C)Wij, 

3 

and if k = both * = and $ = 0, thus F r = and = 0. 

Combining (13) & (14) we get a quartic in the frequency u, one root of which 
is to = 0. This leaves us with: 

u 3 - ig<&u 2 + (2ig^ - l)u - 3ig$ = 0. (15) 

In the non-viscous case (i.e. where ^ = $ = 0) u = 0, ±1, in agreement with 
the solution of (7) and (8). 

If the viscosity (which is proportional to g) is small enough, we expect fre- 
quencies close to the unperturbed values. We therefore set u = 7, 1 + 7 and 
— 1 + 7 where I7I <C 1 and the 7 in each case will be different. 

id = 1 + 7 

Substituting u = 1 + 7 into equation (15) and keeping only the linear terms 
(since I7I <C 1) we find: 

7 = ig(2<$> - ) , (16) 

and 

exp{iiot} = exp{z(l + 7)^} 

= expjzt - gt(2<f> - *)} . 

For stability we require 2$ — ^ > 0. Substituting for $ and ^ we find that 
for all values of k, 2$ — ^ > (see Murray 1994 for this and the following 
two cases). Thus the u = 1 + 7 mode is stable. 
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to = 7 — 1 



In this case we get 

7 = z#(2$ + *) (17) 

and so exp{iiot} shows that for stability we require 2$ + ^ > 0, which holds 
for all values of k. Hence u = 7 — 1 is also stable. 

id = 7 

Here we get: 

7 = -3ig$ . (18) 

This gives exp{iiot} = exp{3g<&t} whose exponent is always positive, and 
hence the u = 7 is an unstable mode for all values of k. Note that in the 
inviscid solution, the u = case is associated with an azimuthal shift with 
radius held constant. 

For large wavelengths (large compared to the particle spacing, which is of 
order h, and thus large wavelength corresponds to hk <C 1), we substitute a 
Gaussian for the kernel and approximate the summation by an integral and 
find: 

Ph , , 

(19) 



/7rA0 

where A0 is the azimuthal particle spacing. Thus 

to = 7 = —3ig<& 



^<«> 2 - < 2 °> 



Since we are assuming that I7I <C 1, we must have: 

ac 

<1, 



to be consistent with the analysis. In the models that we ran, particles were 
equally spaced a scaled distance A0 27r/200 apart, with ac <^ 1 (see 
below for details). Hence, we have found that there is an unstable mode 
when u; <C 1, ie when etc <C 1. 

The time scale of the viscous instability for long wavelengths is: 

Tv « • 

bac 

In the models we ran, a = 0.07, c = 0.05 and e 4/9. Thus 77 1.178 rota- 
tion periods. However, from (20) we see that the most rapid growth occurs 
for short wavelengths for which the approximation (19) to the summation 
is invalid. Hence we expect instabilities with wavelength approximately a 
particle spacing to grow fastest. 
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Figure 2: The onset of viscous instability in a single ring system. 
The figures are at 0.2,4.1 and 5.8 rotational periods. 

3 The Simulations 

We ran an SPH simulation of a single viscous ring of approximately 200 
equally spaced particles orbiting a central massive body. The ring was placed 
at dimensionless radius of 0.9, at which position the ring has a rotational 
period of 5.4 dimensionless time units. As mentioned above, a = 0.07 and 
c = 0.05. Particles started to oscillate very slightly almost immediately and 
after approximately four rotational periods the ring had broadened. The 
particles moved on eccentric orbits and after five rotational periods the ring 
had broadened significantly (see figure 2). It was seen from the simulation 
that the viscous instability sets in after approximately four rotation periods. 

A double ring system was also modelled. If the separation of the rings was 
less than 2h (the radius of influence of the SPH kernel), the two rings acted 
as one and were disturbed in a similar manner as was the one ring system. 
However, if the separation of the rings was greater than 2h, then as usual 
viscosity acts to spread the rings and particles moved onto eccentric orbits. 
The inner rings broke up first after four rotation periods, followed by the outer 
ring at its equivalent rotational phase. The two rings were indistinguishable 
after eight orbits. 



4 Conclusions 

The results of the disk simulations show that SPH gives results in good 
agreement with the standard viscous theory. Deviations from theory eventu- 
ally occur at the inner edge where we observe particles moving from circular 
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orbits. The analysis of Lyubarskij et al. (1994) shows that axisymmetric 
viscous disks can become unstable whenever the viscous drag on fluid in an 
elliptical orbit is a maximum at apastron. This is the case for the standard 
a-disk and it is the case for particles in our simulation. For example, the 
particles in the innermost ring appear to move in and then out which shows 
that they are individually moving on elliptical orbits with apastron near their 
neighbouring ring. Since the viscous force is felt when they move near the 
particles in the neighbouring ring they will be forced into more eccentric 
orbits. 

When modelling viscous disks we found that the innermost ring separated 
from the remaining disk and then became unstable. We then investigated 
the dynamics of a single ring. However, the formulation of SPH is based on 
smoothing. A single ring is the antithesis of a smoothed configuration. The 
single ring smooths itself by rapidly spreading under the viscous forces. The 
combination of analysis and simulation confirms the consistency of SPH. 

However, the resulting effects of the SPH viscous force prescription in ac- 
cretion disk simulations is unphysical. The main problem is that we are 
trying to model a continuum (i.e. a gaseous viscous disk) by a set of discrete 
rings. For disks made up of concentric rings, this problem can be overcome 
if the number of rings used is increased so that the rings will not separate 
due to viscosity by too great a distance (generally less than 2h). As usual 
with particle methods, a compromise must be reached between the number 
of particles and the computation time. Alternatively, one could incorpo- 
rate a variable smoothing length, so that when the distance between rings is 
greater than 2h, the smoothing length increases so that neighbouring rings 
can still communicate. There are still some problems to be overcome with the 
implementation of spatially varying smoothing lengths when sharp density 
gradients are present (see Nelson & Papaloizou 1994). 
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A Appendix 



We want to solve lor 

-S v ? 'j • r ? 'j 

= K'-v.Hr.-r-,). (Al) 

Noting that to first order 

r; = (1 + qi)r io + c;(z X r io ) 

and differentiating with respect to t, (noting dri /dt = z X r 80 , with = 1), 
we find 

v; = (q t - cri)r w + (1 + cr % + ^)(z X r io ), (A2) 

where A denotes a unit vector, r 80 is the unperturbed position vector ol particle 
z, and z is perpendicular to the plane ol the orbit. We can write the terms 
contributing to (Al) as, lor example, 

v; • r; = (qt - <7;)(1 + qi) + ai(&i + 1 + qi) 

^ qi 

Vj ■ r 8 = [(qj - + q^ + (1 + q 3 + <7?)<7;] cos ( 

— (qj + c; - cr 3 ) cos ( - (1 + q t + ?j + <Jj) sin (. (A3) 

where ( = Qj — Qi and r 80 • r J0 = cos(. Similar expressions can be written 
down lor Vj • r 3 and v 4 - • r r Combining these expressions we find 

s = (qi + - cos C) + (&j - °i) sm C- (A4) 

Since < ( <C 1 we can neglect the first term in (A4) to obtain 

s ~ (<Tj — &i) sin (. 



10 



References 

Flebbe,0., Munzel,S., Herold,H., RifFert,H. & Ruder,H. 1994, ApJ, 431, 754 
Lubow,S.H. 1991, ApJ, 381, 268 

Lynden-Bell,D. & Pringle,J.E. 1974, MNRAS, 168, 603 

Lyubarskij,Y.E, Postnov,K.A. & Prokhorov,M.E. 1994, MNRAS, 266, 583 

Monaghan,J.J. 1992, ARA&A, 30, 543 

Murray,J.R. 1994, PhD Thesis, Monash University 

Nelson,R.P. & Papaloizou,J.C.B. 1994, MNRAS, 270, 1 

Pringle,J.E. 1981, ARA&A, 19, 137 



11 



